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Abstract. Modern investigations of magnetic chemically peculiar stars reveal 
a variety of complex processes in their atmospheres. Realistic spectrum synthesis 
modelling of these objects has to take into account anomalous chemical compo- 
sition, strong magnetic field and chemical stratification. These effects complicate 
calculation of theoretical spectra, especially when one has to deal with large num- 
^ ■ ber of lines and wide spectral regions. To overcome the formidable problem of 

comparing model and observed spectra for magnetic chemically peculiar stars, 
a new suite of spectrum synthesis programs was developed. Here we describe in 
detail the synthesis codes, SYNTH3 and SYNTHMAG and present examples of their 
application to various aspects of the peculiar-star surface phenomena. The new 
codes proved to be reliable tools for the line identification, magnetic field deter- 
mination, chemical abundance and stratification analysis. 
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1. Introduction 

Recent major improvements in the quality of spectroscopic observations stimulated a surge of in- 
terest in detailed model atmosphere and chemical abundance studies of magnetic stars. Modern 
spectrographs, such as UVES at the ESO 8-m VLT at Paranal, HARPS at the ESO 3.6-m telescope 
at La Silla, or the NES instrument at the SAO 6-m telescope, can now provide S/N > 300, wide 
spectral coverage echelle spectra of moderately faint stars at the resolving power of A/AA = 4 x 10 4 - 
10 5 . Publicly accessible data archives assist in distributing these spectra to an astronomical com- 
munity far wider than the co-investigators on the original observing proposals or scientists from the 
countries that own the telescopes and instruments. 

Availability of the new observational material has to be matched by the corresponding develop- 
ment of the new analysis techniques, capable of handling large spectral regions and large number of 
lines without compromising accuracy. The complex chemistry of chemically peculiar (CP) stars, cou- 
pled with the prominence of the magnetic field and chemical stratification effects, implies that any 
realistic analysis of the atmospheres of these stars must be partly or entirely based on the method of 
spectrum synthesis. Moreover, to deduce accurate chemical abundances, theoretical computations 
must properly account for the transfer of polarized radiation through the chemically inhomoge- 
neous stellar atmosphere. Numerical schemes implemented in the codes have to be stable against 
rapid variation of the equilibrium model structure and opacity, associated with the chemical clouds, 
located at different atmospheric heights. 

Several computer codes for modelling the time series of Stokes parameter spectra were developed 
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in the past (Wade et al. 2001), but they are typically constrained to the analysis of short wavelength 
regions and are not thoroughly tested with the chemically stratified atmospheres. The lack of general- 
purpose spectrum synthesis codes, suitable for modelling magnetic, chemically stratified stellar 
atmospheres, stimulated development of new tools and major improvement of the existing, widely 
used software. In this paper we describe the current status and recent updates of our two spectrum 
synthesis codes: SYNTHMAG (magnetic spectrum synthesis) and SYNTH3 (computation of non-magnetic 
stellar spectra). The former code represents a major revision and improvement of the early magnetic 
spectrum synthesis code by Piskunov (1999), whereas the latter supersedes the well-known SYNTH 
code (Piskunov 1992). 

2. Physical foundations 

SYNTHMAG and SYNTH3 are designed to calculate spectra emerging from the static, plane-parallel, 
one-dimensional model atmosphere. The Local Thermodynamical Equilibrium assumption is used 
throughout. These approximations are adequate for a wide range of stars in the Main Sequence 
band. In practice, the codes were successfully used to model stars in the range of spectral classes 
from early-B to late-M. 

2.1. Opacities 

The polarized radiation is fully characterized by the Stokes vector I = {/, Q,U, V}^, where the 
Stokes parameters are defined according to Shurcliff (1962). The SYNTHMAG code solves the polarized 
radiative transfer equation in the form 

^— W + J. (1) 

where m denotes the column mass scale, which is used for tabulation of model atmospheres, and 
H = cos 9 is the cosine of the angle between line of sight and the local surface normal. For the sake 
of brevity here we omit complete expressions of the elements of the absorption matrix K, and of the 
emission vector J, since this was discussed in detail by Piskunov & Kochukhov (2002). 

In calculating the metal line opacities, SYNTHMAG takes into account Zeeman or Paschen-Back 
splitting patterns specified in the input line list. The absorption and anomalous dispersion profiles 
for each Zeeman component are described by the Voigt and Faraday- Voigt function, respectively. 
Radiative, van der Waals and Stark broadening is calculated according to the damping parameters 
provided in the line list or using classical expressions (Gray 1992). The chemical and ionization 
equilibrium is calculated for a prescribed vertical distribution of abundances. Absorption profiles 
of all lines contributing to a given wavelength are summed and multiplied by the trigonometric 
functions of the angles 7 and %■> which define orientation of the local field vector with respect to the 
line of sight. The resultant coefficients enter the absorption matrix /C. 

The hydrogen line opacity is determined with the routines developed by Barklem et al. (2000). 
The Stark broadening tables of Stehle (1994) are used to account for the linear Stark effect. Self- 
broadening of hydrogen lines (Barklem et al. 2000) is also included. An optional module permits 
improved calculations of the overlapping lines of the higher hydrogen series members with the 
occupation probability formalism of Daeppen et al. (1987) and Hubeny et al. (1994). 

The Stark broadening parameters of the neutral He lines are computed following Dimitrijevic & 
Sahal-Brechot (1984) and Freudenstein & Cooper (1978). 

The continuous opacity calculation in our spectrum synthesis codes is based upon the routines 
adapted from the ATLAS9 code (Kurucz 1993). We include absorption due to the bound- free and free- 
free transitions of H 1, H^~, H~, He 1, He 11, different metal ions, He - free-free, Rayleigh scattering 
for H 1 and He 1, and electron scattering. 
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Opacity calculations in SYNTH3 are similar to those of SYNTHMAG, except that a non-magnetic 
version of Eq. (TjQ) is used and unpolarized radiation is computed using the usual scalar opacity 
coefficients. It was verified that the line profiles of SYNTH3 are identical to calculations with SYNTHMAG 
if magnetic field strength is set to zero. 

2.2. Ionization and molecular equilibrium 

Ionization and molecular equilibrium in SYNTH3 and SYNTHMAG is obtained with the advanced equa- 
tion of state solver written by N. Piskunov. This code represents a modified and expanded version 
of the routine described in Valenti et al. (1998). The solver includes the treatment of ~ 200 diatomic 
and polyatomic molecules and has been thoroughly tested against similar routines of the MARCS and 
PHOENIX model atmosphere codes for temperatures from ~ 10 5 K down to a few hundred K. 

Atomic partition functions are calculated using the tables given in the PFSAHA subroutine of 
the ATLAS9 code (Kurucz 1993). For a number of heavy elements, especially rare-earths, updated 
partition functions of Cowley & Barisciano (1994) and Sneden (2002) are incorporated. For calcu- 
lations of the ionic populations we typically take into account 6 ionization stages for light elements, 
4 stages for the iron-peak group and 3 stages for heavy elements. 

For diatomic molecules the partition function information is taken from Sauval & Tatum (1984) 
and Irwin (1987). For polyatomic molecules, the data from Irwin (1988) are used. 

2.3. Magnetic field 

Detailed structure of the stellar magnetic field geometry is best studied with the time-resolved 
spectropolarimetric observations using Doppler imaging (Piskunov & Kochukhov 2002) or fitting 
parameters of the multipolar field models (Khalack & Wade 2006). On the other hand, the main 
purpose of SYNTHMAG is to properly account for the Zeeman broadening and splitting in a wide 
wavelength region analysis of the unpolarized spectra of magnetic stars. Information content of a 
few unpolarized stellar observations, available for the majority of CP stars, is usually insufficient to 
constrain even the most simple model of the field geometry, especially when no rotational modulation 
can be detected. This is why the common approach of expanding the field structure in spherical 
harmonics is poorly constrained and one has to resort to a simplified field model, characterized 
by a small number of free parameters. Following this ideology, the SYNTHMAG code calculates the 
local Stokes IQUV spectra for a set of angles between the line of sight and local surface normal 
assuming that i) the field is axisymmetric with respect to the line of sight and ii) the field structure 
is homogeneous and is defined by the three vector components: radial, meridional and azimuthal 
field (see Fig. [1]) . This is essentially equivalent to approximating the field strength distribution with 
a single value - an approach which works surprisingly well for many slowly rotating CP stars. 

3. Numerical details 

Both the SYNTHMAG and SYNTH3 programs are written in Fortran 77, which ensures portability to all 
major computer platforms and straightforward interfacing with alternative software modules. The 
codes use standard BLAS and LAPACK routines, but do not depend on commercial or proprietary 
software libraries. Up to now, the codes were compiled and successfully used on the SunOS and 
HP-UX versions of the Unix operating system, on Mac OSX and Linux. 

In the standard mode, SYNTH3 computes intensity by solving numerically the scalar radiative 
transfer equation for a set of (typically) seven angles. SYNTHMAG proceeds similarly, but solves the 
vector polarized RT equation and derives Stokes IQUV profiles for a given local magnetic field. 

Each code does the following set of operations: 
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Figure 1: Left panel: magnetic field geometry adopted in SYNTHMAG. The field structure is symmetric 
with respect to the line of sight (z-axis). The three magnetic vector components, B r , B a , B m , 
describe the radial, meridional and azimuthal field in the stellar reference frame. Transformation 
of the local field vector to the observer's coordinate system xyz depends on the angle 9 of each of 
the seven annular zones employed for calculation of the disk-integrated profiles. Right panel: typical 
four Stokes parameter disk-center SYNTHMAG calculations for B = 10 kG and T e g = 7700 K. 



• Reading the input model atmosphere (krz format) and atomic/molecular line list in the 
VALD format (Kupka et al. 1999). For SYNTHMAG the input line list also contains information on 
the strength and splitting of the Zeeman components for each line. Depth dependence of chemical 
abundance, magnetic field strength or microturbulent velocity can be specified in the input model 
atmosphere file. 

• Concentrations of ions and molecules are computed with the equation of state module de- 
scribed above. For each layer in the model atmosphere we calculate continuous and line center 
opacities and the Voigt function parameters. The depth grid is refined if necessary. 

• Intensity at several limb angles is computed with the help of formal polarized (or scalar for 
SYNTH3) radiative transfer algorithm, using precomputed line-center opacities. The wavelength grid 
is refined to obtain an accurate description of the line profile shapes. 

• Intensity profiles are broadened for a given projected rotational velocity and macroturbulence. 
The resultant profiles are combined to form the disk-integrated spectra. 

3.1. Integration of the radiative transfer equation 

The polarized radiative transfer equation (1) is solved on discrete vertical grid using the Diagonal 
Element Lambda Operator (DELO) formal solver, introduced by Rees et al. (1989). In our spectrum 
synthesis codes we employ a modified DELO algorithm, which utilizes a parabolic approximation 
of the source function (Socas-Navarro et al. 2000). Piskunov & Kochukhov (2002) have compared 
several widely used polarized radiative transfer algorithms and found that the quadratic DELO 
method is superior in terms of speed and accuracy to the original DELO solver and to all other 
polarized radiative transfer schemes. For typical ATLAS9 model atmospheres of the Main Sequence 
B-F stars the accuracy of Stokes parameter calculation with the DELO algorithm is ~ 10~ 4 . 



5 




_6 5 4 -3 2 1 1 2 6 5-4 3 2 1 1 2 

fc) lj p , 




o.sL_, , , , i , , , , i , , , , i , , , , i , , , , i , , , ,_d 

6147.0 6147.5 6148.0 6148.S 6149.0 6149.5 6150.0 

"Wavelength {A) 

Figure 2: Illustration of the adaptive grid refinement in SYNTHMAG. a) Extra layers are added in 
the model atmosphere if vertical variation of the total opacity is too fast. Large symbols show the 
original discretization of the model; small symbols illustrate the final depth grid, b) Refinement of 
the frequency grid. Initially spectra are calculated for the centers of Zeeman components and for 
the wavelengths in between (large symbols). More points (small symbols) are added until the line 
profile shapes are described with sufficient accuracy. 

Rapid variation of abundances with height in the atmospheres of chemically peculiar stars poses 
further challenges to a radiative transfer solver. The latter has to provide accurate spectra for 
minimum number of layers in the model atmospheres and, at the same time, be robust against steep 
variation of the line opacity resulting from chemical stratification. In the recent paper Kochukhov 
et al. (2006) investigated performance of different scalar radiative transfer algorithms for chemically 
stratified stellar atmospheres. These numerical experiments showed that direct intensity integration 
implemented in Kurucz (1993) codes is remarkably stable and performs well for atmospheres with 
chemical stratification. However, this algorithm cannot be extended to polarized radiation and for 
this reason it is not used in our codes. The quadratic DELO scheme achieves acceptable accuracy of 
~ 10~ 3 for the atmospheres with large abundance gradients if extra vertical grid points are added 
between the layers where opacity changes rapidly. 

3.2. Spectrum synthesis with adaptive grids 

To achieve required precision with the least amount of computations, adaptive grids are extensively 
used in SYNTHMAG and SYNTH3. At the stage of precalculation of the line-center opacities, both 
codes consider the vertical gradient of the monochromatic optical depth. If variation of the optical 
depth exceeds a given threshold, an extra layer is added in the model structure. Tabulated model 
atmosphere quantities are obtained for additional layers with the help of spline interpolation. The 
procedure of the depth grid refinement is essential to obtain accurate spectra for the chemically 
stratified atmospheres. An example of the depth grid refinement for cool Ap-star atmosphere with 
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a strong vertical gradient of the iron abundance is presented in Fig. 2a. 

Most spectrum synthesis codes calculate line profiles on a wavelength grid defined by a fixed 
wavelength step. This leads to computation of excessive number of continuum points with no useful 
information. On the other hand, a coarse equidistant grid may not allow one to resolve details of the 
line profile shapes. SYNTHMAG and SYNTH3 avoid this difficult compromise between the line profile 
accuracy and amount of computations by using adaptive wavelength grid. At the first step, intensity 
is calculated at the centers of spectral lines or Zeeman components, if the magnetic field is present. 
Then spectral points between each pair of the original wavelengths are synthesized. Computed 
intensity is compared with the results of linear interpolation and further frequency points are added 
if the discrepancy between interpolated and computed spectra exceeds a given threshold. In this way, 
more points are calculated for the regions with rapid variation of opacity, whereas the wavelength 
step remains relatively large for continuum regions (see Fig 2b). Using this procedure, we obtain 
theoretical spectra at infinite resolution, performing a factor of 10 to 50 less computations to achieve 
the same accuracy as in the synthesis with a sufficiently fine equidistant wavelength grid. 

3.3. Disk integration 

The final step of the spectrum synthesis calculation with SYNTHMAG and SYNTH3 is integration of the 
local Stokes parameter /intensity spectra to obtain stellar flux profiles. We use an external procedure, 
compatible with both the SYNTHMAG and SYNTH3 output files, to perform the disk integration. 

In the default mode, the local intensity and Stokes parameter spectra are calculated for seven 
annular surface zones, chosen in such a way that their projected surface areas are equal (see left 
panel in Fig. 1). Numerical experiments have demonstrated that this discretization of the stellar 
surface is sufficient for computation of high-precision unpolarized spectra. On the other hand, due 
to substantial variation of the shape and intensity of the Stokes QUV local profiles over the stellar 
surface in any realistic magnetic field configuration, co-addition of seven linear and circular polariza- 
tion local spectra provides only a rough estimate of the signal expected in the QUV disk-integrated 
spectra. 

Intensity spectra are convolved with the appropriate rotational and radial-tangential macrotur- 
bulence broadening profiles (Gray 1992) and then added with equal weights. This assumes spherically 
symmetric star, covered by a homogeneous magnetic field. We note, however, that SYNTHMAG and 
SYNTH3 are, in fact, independent of the latter assumptions of the disk integration and can be used to 
supply intensity spectra for a more sophisticated disk integration procedure. For instance, different 
sets of intensity spectra can be combined to simulate observations of rotationally distorted stars, 
stars covered with chemical, temperature or magnetic spots, etc. Below we give an example of such 
an advanced usage of SYNTHMAG for calculation of the Zeeman resolved spectral lines corresponding 
to dipolar field geometry. 

3.4. Reconstruction of chemical stratification 

Early attempts to deduce vertical chemical stratification in the atmospheres of Ap stars were based 
on the trial-and-error method (Ryabchikova et al. 2002). Abundances in the upper and lower parts 
of the stellar atmosphere were manually adjusted to fit the line profile shapes and to remove a 
discrepancy between the strength of lines of different excitation potential, equivalent width and 
ionization stages. We have recently developed an automatic procedure, DDAFit, to find chemical 
abundance gradients from observed stellar spectra (Ryabchikova et al. 2005). The DDAFit script is 
written in IDL and provides an optimization and visualization interface to the spectrum synthesis 
calculations with SYNTHMAG and SYNTH3. Vertical abundance distributions are described with the four 
parameters: chemical abundance in the upper atmosphere, abundance in deep layers, the vertical 
position of abundance step and the width of the transition region where chemical abundance changes 
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Figure 3: Examples of the SYNTHMAG modelling of the spectra of magnetic CP stars. Top panel: 
comparison between observations of the red Si n doublet at A 6347 and 6371 A and spectrum 
synthesis for HD 137509 - the star with second-largest field known ((B) = 29 kG, Kochukhov 
2006). Bottom left panel: detection of the Zeeman resolved components in the Fe n 5018 A line in 
the spectrum of cool Ap star HD 178892 (Ryabchikova et al. 2006a). In this plot solid lines show 
observations obtained with the NES instrument at the SAO 6-m telescope, whereas dashed line 
shows SYNTHMAG calculation for (B) = 17.5 kG. Bottom right panel: observations (histogram) of 
the Y II 5087 A line in the strongly magnetic cool Ap star HD 154708. Usual SYNTHMAG calculation 
for (B) = 24 kG (dashed line) is compared with the result of combining spectra for a dipolar-like 
distribution of the field strength in the range from 29 to 16.5 kG (solid line). 



between the two values. All four parameters can be optimized simultaneously with the Levenberg- 
Marquardt least-squares fitting routine (Bevington & Robinson 1992) and based on observations of 
unlimited number of spectral regions, possibly using different weights in accordance to the quality 
or relative importance of the observations of particular spectral features. The program derives one 
chemical stratification profile at a time, but is able to account for any number of fixed stratified 
abundances. 
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4. Applications 

4.1. Magnetic stars 

SYMTHMAG was used in a number of recent studies of magnetic Ap stars. In such analyses polarized 
radiative transfer calculations are essential to derive accurate abundances, resolve contributors to 
complex blends and to measure magnetic field strength. 

Fig. 3 illustrates application of SYNTHMAG to three different chemically peculiar stars with strong 
magnetic field. For the Si-peculiar Bp star HD 137509 (Kochukhov 2006) we show comparison of 
observations and theoretical spectra in the region of the Si II doublet at A 6347 and 6371 A. Strikingly 
different profile shapes of the two, otherwise similar, lines can be traced to the different Zeeman 
splitting structure. Resolved Zeeman components in the Si n 6371 A line indicate (B) = 29 kG, 
which is a second-largest field ever observed in a non-degenerate star. This very large field strength 
is confirmed by the spectrum synthesis of other metal lines. 

We also show application of SYNTHMAG to the analysis of HD 178892 - a star discovered by 
the magnetic surveys at the SAO 6-m telescope (Elkin et al. 2002). The follow up high-resolution 
spectroscopy with NES (Ryabchikova et al. 2006a) reveal resolved Zeeman components in many 
lines, for example Fe n 5018 A shown in Fig. 3. Low T c g, strong magnetic field ((B) = 17.5 kG) 
and relatively short rotation period make HD 178892 an excellent candidate for a detailed analysis 
of the horizontal structure of magnetic field and chemical abundance spots. 

HD 154708 is another unique cool Ap star, showing very strong magnetic field ((B) = 24.5 kG, 
Hubrig et al. 2005). Recent time-series observations at ESO suggest the presence of low-amplitude 
p-mode pulsations in this star (Kurtz et al. 2006). Due to an extremely slow rotation, Zeeman 
components are fully resolved for many spectral lines in HD 154708. The field is so strong however, 
that the standard SYNTHMAG approximation of the surface field strength distribution with a single 
value of the field modulus fails (see Fig. 3). This is evident from the comparison of spectrum 
synthesis and observed profiles of the a components in the Zeeman split triplet lines: a components 
are broader and shallower in observations, indicating a non-negligible surface scatter of the field 
strength. This situation can be addressed with a straightforward extension of the spectrum synthesis 
with SYNTHMAG. Fig. 3 shows that a much better agreement with observations can be achieved if the 
stellar flux profiles are produced from a set of SYNTHMAG computations for different field strengths, 
corresponding to the B p = 30 kG dipolar field topology. 

4.2. Chemically stratified atmospheres 

Both SYNTHMAG and SYNTH3 can calculate spectra for the prescribed inhomogeneous vertical dis- 
tribution of chemical elements. Parameters of the simple step-function distribution of abundances 
with depth can be adjusted with the help of the DDAFit script described above. 

A number of recent chemical abundance analyses of Ap stars accounted for chemical stratification 
in the spectrum synthesis with SYNTHMAG and SYNTH3 and obtained vertical distributions of selected 
elements using the DDAFit procedure. Fig. 4 illustrates application of SYNTHMAG to the problem of 
deriving Fe stratification in the strongly magnetic Ap star HD 144897 (Ryabchikova et al. 2006b). 
Another example of employing SYNTH3 and DDAFit to the chemical stratification modelling can be 
found in Ryabchikova et al. (2005). 

5. Conclusions 

We have developed a new family of spectrum synthesis codes for modelling the spectra of stars 
with various atmospheric parameters and chemistry. The codes SYNTHMAG and SYNTH3 can be used 
to calculate line profiles of normal stars with different T e g, as well as peculiar stellar objects with 
strong magnetic fields, non-solar chemical composition and vertical chemical abundance gradients. 
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Figure 4: Spectrum synthesis calculation for chemically stratified atmosphere of the Ap star 
HD 144897 (Ryabchikova et al. 2006b). Symbols show observed line profiles. Dashed curve illus- 
trates theoretical spectrum for constant Fe abundance. A better fit to observations (solid line) 
could be achieved by taking into account vertical stratification of Fe. 

Inclusion of the molecular opacities and employment of the general equation of state solver permits 
application of our codes to a wide variety of stars across the H-R diagramm. 
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